Isomap — Nonlinear Dimensionality Reduction (Geodesic Distances)#
Isomap (Isometric Mapping) is a manifold learning method that tries to preserve walking distances along a surface, not straight-line “as the crow flies” distances.
What you’ll learn#
why Euclidean distance can be misleading on curved manifolds (Swiss roll)
how Isomap works: kNN graph → geodesic distances → classical MDS
how to implement Isomap from scratch (NumPy + shortest paths)
how to diagnose common failure modes (noise, disconnected graphs)
import heapq
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
from scipy.linalg import orthogonal_procrustes
from sklearn.manifold import Isomap as SklearnIsomap
from sklearn.metrics import pairwise_distances
from sklearn.neighbors import NearestNeighbors
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)
Prerequisites#
Euclidean distance + basic linear algebra
Graphs and shortest paths (Dijkstra / Floyd–Warshall)
Eigendecomposition (classical MDS)
1) Intuition: “Measuring distance by walking paths, not flying”#
If you’re in a city with rivers and bridges, two places can be:
close “as the crow flies” (straight line)
but far if you must walk along streets (path distance)
Isomap is built on the same idea:
Euclidean distance ≈ “flying”
geodesic distance ≈ “walking along the manifold”
Swiss-roll metaphor#
A Swiss roll is a 2D sheet rolled into 3D.
Points on different layers can be close in 3D Euclidean distance.
But along the sheet, they can be far.
Isomap tries to “unroll” the sheet by preserving geodesic distances.
def make_swiss_roll(n_samples: int = 350, noise: float = 0.0, random_state: int = 7):
rng_local = np.random.default_rng(random_state)
# Roughly matches sklearn's swiss-roll parameterization
t = 1.5 * np.pi * (1.0 + 2.0 * rng_local.random(n_samples)) # [1.5π, 4.5π]
h = 21.0 * rng_local.random(n_samples)
x = t * np.cos(t)
z = t * np.sin(t)
y = h
X = np.column_stack([x, y, z])
if noise > 0:
X = X + rng_local.normal(scale=noise, size=X.shape)
return X, t, h
X, t, h = make_swiss_roll(n_samples=350, noise=0.0, random_state=7)
X.shape
(350, 3)
fig = px.scatter_3d(
x=X[:, 0],
y=X[:, 1],
z=X[:, 2],
color=t,
title="Swiss roll in 3D (color = intrinsic coordinate t)",
labels={"x": "x", "y": "y (height)", "z": "z", "color": "t"},
)
fig.update_traces(marker=dict(size=4, opacity=0.85))
fig.update_layout(scene=dict(aspectmode="data"))
fig.show()
fig = px.scatter(
x=t,
y=h,
color=t,
title="Ground truth 2D coordinates (t, height)",
labels={"x": "t (along the sheet)", "y": "height", "color": "t"},
)
fig.update_traces(marker=dict(size=6, opacity=0.85))
fig.show()
2) Graph-based explanation#
Isomap is a graph + geometry method:
Build a k-nearest neighbors (kNN) graph in the original space.
Use shortest paths in that graph to approximate geodesic distances.
Run classical MDS on those geodesic distances to get a low-dimensional embedding.
2.1 kNN graph#
each point is a node
connect each node to its k nearest neighbors (by Euclidean distance)
edge weight = that local Euclidean distance
Locally, Euclidean distance is a good approximation to distance along the manifold.
def knn_graph(X: np.ndarray, k: int) -> list[list[tuple[int, float]]]:
"""Undirected kNN graph with edge weights = Euclidean distances."""
nn = NearestNeighbors(n_neighbors=k + 1, algorithm="auto")
nn.fit(X)
dists, idx = nn.kneighbors(X)
idx = idx[:, 1:]
dists = dists[:, 1:]
n = X.shape[0]
nbr_dicts: list[dict[int, float]] = [dict() for _ in range(n)]
for i in range(n):
for j, d in zip(idx[i], dists[i]):
j = int(j)
w = float(d)
nbr_dicts[i][j] = min(nbr_dicts[i].get(j, np.inf), w)
nbr_dicts[j][i] = min(nbr_dicts[j].get(i, np.inf), w)
return [sorted(nbrs.items(), key=lambda t: t[1]) for nbrs in nbr_dicts]
def connected_components(neighbors: list[list[tuple[int, float]]]):
n = len(neighbors)
comp = -np.ones(n, dtype=int)
cid = 0
for start in range(n):
if comp[start] != -1:
continue
stack = [start]
comp[start] = cid
while stack:
u = stack.pop()
for v, _ in neighbors[u]:
if comp[v] == -1:
comp[v] = cid
stack.append(v)
cid += 1
return comp, cid
def ensure_connected_knn_graph(X: np.ndarray, k_start: int = 4, k_max: int = 30):
for k in range(k_start, k_max + 1):
neighbors = knn_graph(X, k=k)
_, n_comp = connected_components(neighbors)
if n_comp == 1:
return neighbors, k
raise ValueError(f"Graph stayed disconnected for k in [{k_start}, {k_max}]")
k_target = 10
neighbors, k = ensure_connected_knn_graph(X, k_start=k_target, k_max=30)
k
10
# Visualize the kNN graph on a smaller subset (edges get cluttered fast)
n_vis = 160
idx_vis = rng.choice(X.shape[0], size=n_vis, replace=False)
Xv = X[idx_vis]
tv = t[idx_vis]
neighbors_v, k_v = ensure_connected_knn_graph(Xv, k_start=10, k_max=25)
xs, ys, zs = [], [], []
for a, nbrs in enumerate(neighbors_v):
for b, _w in nbrs:
if a < b:
xs.extend([Xv[a, 0], Xv[b, 0], None])
ys.extend([Xv[a, 1], Xv[b, 1], None])
zs.extend([Xv[a, 2], Xv[b, 2], None])
fig = go.Figure()
fig.add_trace(
go.Scatter3d(
x=xs,
y=ys,
z=zs,
mode="lines",
line=dict(color="rgba(0,0,0,0.18)", width=1),
name=f"kNN edges (k={k_v})",
hoverinfo="skip",
)
)
fig.add_trace(
go.Scatter3d(
x=Xv[:, 0],
y=Xv[:, 1],
z=Xv[:, 2],
mode="markers",
marker=dict(size=4, color=tv, colorscale="Viridis", opacity=0.9, colorbar=dict(title="t")),
name="points",
)
)
fig.update_layout(
title=f"kNN graph on Swiss roll (n={n_vis}, k={k_v})",
scene=dict(aspectmode="data"),
)
fig.show()
2.2 Geodesic distances via shortest paths#
Once we have a neighborhood graph, we approximate geodesic distances by computing shortest path distances in the graph:
Use Dijkstra from every source node (common when edges are non-negative).
Or use Floyd–Warshall (simple, but \(O(n^3)\); good for small graphs).
These shortest-path distances are our estimate of “walking distance” along the manifold.
def dijkstra_distances(neighbors: list[list[tuple[int, float]]], source: int):
n = len(neighbors)
dist = np.full(n, np.inf)
dist[source] = 0.0
heap: list[tuple[float, int]] = [(0.0, source)]
while heap:
d_u, u = heapq.heappop(heap)
if d_u != dist[u]:
continue
for v, w_uv in neighbors[u]:
nd = d_u + w_uv
if nd < dist[v]:
dist[v] = nd
heapq.heappush(heap, (nd, v))
return dist
def dijkstra_with_prev(neighbors: list[list[tuple[int, float]]], source: int):
n = len(neighbors)
dist = np.full(n, np.inf)
prev = -np.ones(n, dtype=int)
dist[source] = 0.0
heap: list[tuple[float, int]] = [(0.0, source)]
while heap:
d_u, u = heapq.heappop(heap)
if d_u != dist[u]:
continue
for v, w_uv in neighbors[u]:
nd = d_u + w_uv
if nd < dist[v]:
dist[v] = nd
prev[v] = u
heapq.heappush(heap, (nd, v))
return dist, prev
def all_pairs_shortest_paths(neighbors: list[list[tuple[int, float]]]):
n = len(neighbors)
D = np.empty((n, n), dtype=float)
for s in range(n):
D[s] = dijkstra_distances(neighbors, source=s)
return D
def adjacency_matrix(neighbors: list[list[tuple[int, float]]]):
n = len(neighbors)
W = np.full((n, n), np.inf)
np.fill_diagonal(W, 0.0)
for i, nbrs in enumerate(neighbors):
for j, w in nbrs:
if w < W[i, j]:
W[i, j] = w
return W
def floyd_warshall(W: np.ndarray):
D = W.copy()
n = D.shape[0]
for k in range(n):
D = np.minimum(D, D[:, [k]] + D[[k], :])
return D
D_geo = all_pairs_shortest_paths(neighbors)
D_euc = pairwise_distances(X)
np.isfinite(D_geo).all(), D_geo.shape
(True, (350, 350))
# A quick sanity check that Floyd–Warshall matches Dijkstra (on a smaller graph)
idx_small = rng.choice(X.shape[0], size=120, replace=False)
X_small = X[idx_small]
neighbors_small, k_small = ensure_connected_knn_graph(X_small, k_start=10, k_max=25)
D_dij = all_pairs_shortest_paths(neighbors_small)
W_small = adjacency_matrix(neighbors_small)
D_fw = floyd_warshall(W_small)
np.max(np.abs(D_dij - D_fw))
7.105427357601002e-15
3) Plotly visualizations#
3.1 Euclidean vs geodesic distances (matrix view)#
If we sort points by the intrinsic coordinate \(t\), a Swiss roll has a clear “rolled up” structure.
Euclidean distances can connect across layers.
Geodesic distances tend to respect the manifold.
order = np.argsort(t)
m = 120
idx_h = order[:m]
E = D_euc[np.ix_(idx_h, idx_h)]
G = D_geo[np.ix_(idx_h, idx_h)]
fig = make_subplots(
rows=1,
cols=2,
subplot_titles=("Euclidean distance (flying)", "Geodesic distance (walking)"),
)
fig.add_trace(go.Heatmap(z=E, colorscale="Viridis", showscale=False), row=1, col=1)
fig.add_trace(go.Heatmap(z=G, colorscale="Viridis", showscale=True), row=1, col=2)
fig.update_layout(title=f"Euclidean vs geodesic distances (subset m={m}, sorted by t)")
fig.show()
3.2 Two points: “close if you fly”, “far if you walk”#
We’ll automatically find two points that have:
small Euclidean distance, but
large difference in intrinsic coordinate \(t\) (likely different layers)
Then we compare:
the straight-line chord (Euclidean)
the shortest path along the kNN graph (geodesic estimate)
def find_cross_layer_pair(D_euc: np.ndarray, t: np.ndarray, min_delta_t: float = 7.0):
D = D_euc.copy()
np.fill_diagonal(D, np.inf)
for thr in [min_delta_t, 6.0, 5.0]:
mask = np.abs(t[:, None] - t[None, :]) > thr
cand = np.where(mask, D, np.inf)
i, j = np.unravel_index(np.argmin(cand), cand.shape)
if np.isfinite(cand[i, j]):
return int(i), int(j), float(thr)
raise ValueError("No suitable cross-layer pair found")
i, j, thr = find_cross_layer_pair(D_euc, t, min_delta_t=7.0)
dist_i, prev_i = dijkstra_with_prev(neighbors, source=i)
path = [j]
cur = j
while cur != i and cur != -1:
cur = int(prev_i[cur])
path.append(cur)
path = path[::-1]
print(f"Chosen pair: i={i}, j={j} (|Δt| > {thr})")
print(f"Euclidean distance: {D_euc[i, j]:.3f}")
print(f"Geodesic distance (graph shortest path): {D_geo[i, j]:.3f}")
print(f"Number of edges in path: {len(path) - 1}")
Chosen pair: i=176, j=196 (|Δt| > 7.0)
Euclidean distance: 8.904
Geodesic distance (graph shortest path): 25.947
Number of edges in path: 5
path_xyz = X[path]
fig = go.Figure()
fig.add_trace(
go.Scatter3d(
x=X[:, 0],
y=X[:, 1],
z=X[:, 2],
mode="markers",
marker=dict(size=3, color=t, colorscale="Viridis", opacity=0.55),
name="points",
)
)
fig.add_trace(
go.Scatter3d(
x=[X[i, 0], X[j, 0]],
y=[X[i, 1], X[j, 1]],
z=[X[i, 2], X[j, 2]],
mode="lines",
line=dict(color="red", width=7),
name="Euclidean chord",
)
)
fig.add_trace(
go.Scatter3d(
x=path_xyz[:, 0],
y=path_xyz[:, 1],
z=path_xyz[:, 2],
mode="lines+markers",
line=dict(color="black", width=4),
marker=dict(size=4, color="black"),
name="Graph shortest path",
)
)
fig.add_trace(
go.Scatter3d(
x=[X[i, 0]],
y=[X[i, 1]],
z=[X[i, 2]],
mode="markers",
marker=dict(size=9, color="red"),
name="start",
)
)
fig.add_trace(
go.Scatter3d(
x=[X[j, 0]],
y=[X[j, 1]],
z=[X[j, 2]],
mode="markers",
marker=dict(size=9, color="red"),
name="end",
)
)
fig.update_layout(
title=f"Close in Euclidean, far in geodesic (k={k})",
scene=dict(aspectmode="data"),
)
fig.show()
4) Classical MDS on geodesic distances (the “Isomap” embedding step)#
Classical MDS takes a distance matrix \(D\) and finds points in a low-dimensional Euclidean space whose pairwise distances match \(D\) as closely as possible.
Given distances \(D_{ij}\), classical MDS:
squares distances: \(D^2\)
double-centers them to build an inner-product (Gram) matrix \(B\)
takes the top eigenvectors/eigenvalues of \(B\)
Isomap is simply:
classical MDS, but using geodesic distances instead of Euclidean ones.
Implementation note: to match the supervised-learning notebooks, we’ll wrap the full pipeline in a small ScratchIsomap class (fit → geodesic distances → MDS embedding).
def classical_mds(D: np.ndarray, n_components: int = 2):
if not np.isfinite(D).all():
raise ValueError("Distance matrix has non-finite entries (graph disconnected?)")
n = D.shape[0]
D2 = D**2
J = np.eye(n) - np.ones((n, n)) / n
B = -0.5 * J @ D2 @ J
eigvals, eigvecs = np.linalg.eigh(B)
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
eigvecs = eigvecs[:, idx]
eigvals = np.maximum(eigvals, 0.0)
Y = eigvecs[:, :n_components] * np.sqrt(eigvals[:n_components])
return Y, eigvals
class ScratchIsomap:
def __init__(
self,
n_neighbors: int = 10,
n_components: int = 2,
shortest_path: str = "dijkstra",
ensure_connected: bool = True,
k_max: int = 30,
):
self.n_neighbors = int(n_neighbors)
self.n_components = int(n_components)
self.shortest_path = str(shortest_path)
self.ensure_connected = bool(ensure_connected)
self.k_max = int(k_max)
def fit(self, X: np.ndarray):
X = np.asarray(X, dtype=float)
if X.ndim != 2:
raise ValueError('X must be a 2D array')
n = X.shape[0]
if n < 2:
raise ValueError('X must have at least 2 samples')
if self.n_components < 1:
raise ValueError('n_components must be >= 1')
if self.n_neighbors < 1:
raise ValueError('n_neighbors must be >= 1')
if self.n_neighbors >= n:
raise ValueError('n_neighbors must be < n_samples')
if self.ensure_connected:
neighbors, k_used = ensure_connected_knn_graph(X, k_start=self.n_neighbors, k_max=self.k_max)
else:
neighbors = knn_graph(X, k=self.n_neighbors)
k_used = self.n_neighbors
if self.shortest_path == 'dijkstra':
D_geo = all_pairs_shortest_paths(neighbors)
elif self.shortest_path in {'floyd-warshall', 'floyd_warshall'}:
D_geo = floyd_warshall(adjacency_matrix(neighbors))
else:
raise ValueError("shortest_path must be 'dijkstra' or 'floyd_warshall'")
Y, eigvals = classical_mds(D_geo, n_components=self.n_components)
self.X_ = X
self.neighbors_ = neighbors
self.n_neighbors_ = k_used
self.geodesic_distances_ = D_geo
self.embedding_ = Y
self.eigvals_ = eigvals
return self
def fit_transform(self, X: np.ndarray) -> np.ndarray:
return self.fit(X).embedding_
scratch_iso = ScratchIsomap(n_neighbors=k, n_components=2, shortest_path='dijkstra', ensure_connected=False).fit(X)
# Use the scratch model outputs downstream
neighbors = scratch_iso.neighbors_
k = scratch_iso.n_neighbors_
D_geo = scratch_iso.geodesic_distances_
Y = scratch_iso.embedding_
eigvals = scratch_iso.eigvals_
eigvals[:10]
array([42278.283 , 38156.4554, 23052.3745, 16412.0213, 4571.1708,
4122.0284, 2861.4723, 2175.4576, 1825.6325, 1413.0827])
fig = px.line(
x=np.arange(1, 11),
y=eigvals[:10],
markers=True,
title="Top eigenvalues (classical MDS on geodesic distances)",
labels={"x": "component", "y": "eigenvalue"},
)
fig.show()
fig = make_subplots(rows=1, cols=2, subplot_titles=("Ground truth (t, h)", "Isomap embedding (from scratch)"))
fig.add_trace(
go.Scatter(
x=t,
y=h,
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
showlegend=False,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=Y[:, 0],
y=Y[:, 1],
mode="markers",
marker=dict(
size=4,
color=t,
colorscale="Viridis",
opacity=0.8,
colorbar=dict(title="t"),
),
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_yaxes(title_text="h (height)", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)
fig.update_layout(title="Unrolling the Swiss roll: Isomap vs ground truth")
fig.show()
5) Algorithm steps (summary)#
Neighborhood graph
build a kNN graph (nodes = points, weighted edges = local Euclidean distances)
Shortest paths
compute all-pairs shortest paths to approximate geodesic distances
common choices: Dijkstra (run from each node) or Floyd–Warshall (small graphs)
Embedding
run classical MDS on the geodesic distance matrix to get low-dimensional coordinates
5.1 Swiss-roll unfolding animation (Plotly)#
To make the “unrolling” feel concrete, we’ll animate a transition from:
the original 3D Swiss roll, to
a 2D Isomap embedding placed onto a flat plane
This is just a visual interpolation, not part of the algorithm.
def zscore(A: np.ndarray):
return (A - A.mean(axis=0)) / A.std(axis=0)
# Align the Isomap embedding to the "true" unrolled coordinates for a nicer animation.
T_true = np.column_stack([t, h])
Y0 = zscore(Y)
T0 = zscore(T_true)
R, _ = orthogonal_procrustes(Y0, T0)
Y_aligned = Y0 @ R
Xc = X - X.mean(axis=0)
scale = float(np.mean(Xc.std(axis=0)))
plane = np.column_stack([Y_aligned[:, 0], Y_aligned[:, 1], np.zeros(X.shape[0])]) * scale
mins = np.minimum(Xc.min(axis=0), plane.min(axis=0))
maxs = np.maximum(Xc.max(axis=0), plane.max(axis=0))
pad = 0.06 * (maxs - mins)
ranges = [(mins[i] - pad[i], maxs[i] + pad[i]) for i in range(3)]
alphas = np.linspace(0.0, 1.0, 31)
def frame_points(alpha: float):
P = (1.0 - alpha) * Xc + alpha * plane
return P
P0 = frame_points(alphas[0])
trace0 = go.Scatter3d(
x=P0[:, 0],
y=P0[:, 1],
z=P0[:, 2],
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.9, cmin=t.min(), cmax=t.max(), colorbar=dict(title="t")),
)
frames = []
for a in alphas:
P = frame_points(float(a))
frames.append(
go.Frame(
data=[
go.Scatter3d(
x=P[:, 0],
y=P[:, 1],
z=P[:, 2],
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.9, cmin=t.min(), cmax=t.max()),
)
],
name=f"{a:.2f}",
)
)
steps = [
{
"args": [[fr.name], {"frame": {"duration": 0, "redraw": True}, "mode": "immediate"}],
"label": fr.name,
"method": "animate",
}
for fr in frames
]
fig = go.Figure(data=[trace0], frames=frames)
fig.update_layout(
title="Swiss-roll unfolding (interpolating to a 2D Isomap embedding plane)",
scene=dict(
xaxis=dict(range=ranges[0], title="x"),
yaxis=dict(range=ranges[1], title="y"),
zaxis=dict(range=ranges[2], title="z"),
aspectmode="data",
),
updatemenus=[
{
"type": "buttons",
"showactive": False,
"x": 0.05,
"y": 1.08,
"buttons": [
{
"label": "Play",
"method": "animate",
"args": [
None,
{
"frame": {"duration": 80, "redraw": True},
"fromcurrent": True,
"transition": {"duration": 0},
},
],
},
{
"label": "Pause",
"method": "animate",
"args": [[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate"}],
},
],
}
],
sliders=[
{
"active": 0,
"currentvalue": {"prefix": "unfold α="},
"pad": {"t": 45},
"steps": steps,
}
],
)
fig.show()
6) Failure modes#
6.1 Sensitivity to noise#
Isomap relies on the kNN graph approximating the true manifold locally.
With too much noise, local neighborhoods get corrupted.
That makes geodesic distances (shortest paths) unreliable.
6.2 Disconnected graphs#
If k is too small, the kNN graph can become disconnected.
Some geodesic distances become infinite.
Classical MDS can’t embed an infinite distance matrix.
# Noise example
Xn, tn, hn = make_swiss_roll(n_samples=350, noise=1.0, random_state=7)
neighbors_n, k_n = ensure_connected_knn_graph(Xn, k_start=10, k_max=30)
D_geo_n = all_pairs_shortest_paths(neighbors_n)
Yn, _ = classical_mds(D_geo_n, n_components=2)
fig = make_subplots(rows=1, cols=2, subplot_titles=("Isomap (clean)", "Isomap (noise=1.0)"))
fig.add_trace(
go.Scatter(
x=Y[:, 0],
y=Y[:, 1],
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
showlegend=False,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=Yn[:, 0],
y=Yn[:, 1],
mode="markers",
marker=dict(size=4, color=tn, colorscale="Viridis", opacity=0.8, colorbar=dict(title="t")),
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(title_text="component 1", row=1, col=1)
fig.update_yaxes(title_text="component 2", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)
fig.update_layout(title="Noise can distort neighborhoods and geodesic distances")
fig.show()
# Disconnected graph example (k too small)
neighbors_bad = knn_graph(X, k=3)
comp_bad, n_comp_bad = connected_components(neighbors_bad)
sizes = np.bincount(comp_bad)
print(f"Connected components: {n_comp_bad}")
print("Component sizes:", sizes)
fig = px.scatter_3d(
x=X[:, 0],
y=X[:, 1],
z=X[:, 2],
color=comp_bad.astype(str),
title="kNN graph can disconnect when k is too small (here: k=3)",
labels={"x": "x", "y": "y", "z": "z", "color": "component"},
)
fig.update_traces(marker=dict(size=4, opacity=0.85))
fig.update_layout(scene=dict(aspectmode="data"))
fig.show()
D_geo_bad = all_pairs_shortest_paths(neighbors_bad)
finite_fraction = float(np.isfinite(D_geo_bad).mean())
print(f"Finite geodesic distance fraction: {finite_fraction:.3f}")
if not np.isfinite(D_geo_bad).all():
print("Disconnected graph → some geodesic distances are infinite. Increase k (or use a radius graph).")
Connected components: 2
Component sizes: [344 6]
Finite geodesic distance fraction: 0.966
Disconnected graph → some geodesic distances are infinite. Increase k (or use a radius graph).
7) Practical usage: scikit-learn’s Isomap#
scikit-learn’s Isomap does the same conceptual steps (kNN graph → shortest paths → embedding), but with a well-tested implementation.
iso = SklearnIsomap(n_neighbors=k, n_components=2)
Y_sklearn = iso.fit_transform(X)
fig = make_subplots(rows=1, cols=2, subplot_titles=("From scratch", "sklearn"))
fig.add_trace(
go.Scatter(
x=Y[:, 0],
y=Y[:, 1],
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8),
showlegend=False,
),
row=1,
col=1,
)
fig.add_trace(
go.Scatter(
x=Y_sklearn[:, 0],
y=Y_sklearn[:, 1],
mode="markers",
marker=dict(size=4, color=t, colorscale="Viridis", opacity=0.8, colorbar=dict(title="t")),
showlegend=False,
),
row=1,
col=2,
)
fig.update_xaxes(title_text="component 1", row=1, col=1)
fig.update_yaxes(title_text="component 2", row=1, col=1)
fig.update_xaxes(title_text="component 1", row=1, col=2)
fig.update_yaxes(title_text="component 2", row=1, col=2)
fig.update_layout(title=f"Isomap embeddings (k={k})")
fig.show()
Exercises#
Try different values of
k. When does the embedding start to “short-circuit” (connect across layers)?Increase
noisegradually. How much noise can Isomap tolerate before it fails to unroll the Swiss roll?Replace the kNN graph with an ε-radius graph. How does that change connectivity and the embedding?
References#
Tenenbaum, de Silva, Langford (2000): A Global Geometric Framework for Nonlinear Dimensionality Reduction
scikit-learn docs:
sklearn.manifold.Isomap